CO2 valence (orbs 8 - 11) ePS results: MFPADs and pol geom dependence¶

22/08/24

Quick revisit of old ePS results (2019) to compare with ELI measurements.

Working from recent N2O notebooks, but may still need some plotter updates. (See https://phockett.github.io/ePSdata/N2O-preliminary/N2O_orbs6-9_LF-MF_preliminary_220224-tidy.html.)

For pol state/rotations, see https://epsproc.readthedocs.io/en/3d-afpad-dev/special_topics/frame_definitions_docs_130323.html#Pol-with-the-Epol-class and https://phockett.github.io/ePSdata/N2O-preliminary/N2O_orbs6-9_MF-polStates_180324.html.

20/03/24

Added basic MF-PAD results (linear polarization, $(x,y,z)$ alignments).

22/02/24 PH

1st processing for N2O results.

For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

Setup¶

In [1]:
%matplotlib inline

# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

/tmp/ipykernel_80233/843188395.py:7: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML
In [2]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr

import matplotlib.pyplot as plt

from datetime import datetime as dt
timeString = dt.now()

import epsproc as ep

# Plotters
from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob

hvPlotters.setPlotters(width = 700, snsStyle='whitegrid')

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 

* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.

* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 

* Set Holoviews with bokeh.
In [3]:
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns

sns.set(rc={'figure.figsize':(10,6)})  # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid")  # Set plot style
sns.set_palette("Paired")   # Set colour mapping

# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')

/tmp/ipykernel_80233/1412183033.py:14: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`
  set_matplotlib_formats('svg', 'pdf')
In [4]:
# xr.set_options(display_style='html')

Load data¶

In [5]:
# import warnings
# # warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
# warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [6]:
# # Scan for subdirs, based on existing routine in getFiles()

# fileBase = Path('/home/paul/ePS/OCS/OCS_survey')  # Data dir on Stimpy

# fileBase = Path('/home/paul/fock-mount/globalhome/eps/N2O/N2O_valence') # Data dir on Jake

fileBase = Path('/home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2') # Data dir on Jake via Docker QM3 container
In [7]:
# TODO: fix orb label here, currently relies on (different) fixed format

data = ePSmultiJob(fileBase, verbose = 0, ext='.inp.out', jobStructure='subDirs')  # TODO: fix subDirs usage!
                                                                                    # CO2 case: fix by setting .inp.out, otherwise also reads base dir (older results)

data.scanFiles(outputKeyType='orb')
data.jobsSummary()

Found 4 directories, with 4 files.

*** Job orb10 details
Key: orb10
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_1pg_0.1-50.1eV, 1 file(s).
{   'batch': 'ePS co2, batch co2_1pg_0.1-50.1eV, orbital E1G',
    'event': ' CO2 1pig-1',
    'orbE': -14.871022583432442,
    'orbLabel': 'orb10',
    'symmetryLabel': '# CO2 1pig-1'}

*** Job orb8 details
Key: orb8
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_1pu_0.1-50.1eV, 1 file(s).
{   'batch': 'ePS co2, batch co2_1pu_0.1-50.1eV, orbital E1U',
    'event': ' CO2 1piu-1',
    'orbE': -19.79084121670707,
    'orbLabel': 'orb8',
    'symmetryLabel': '# CO2 1piu-1'}

*** Job orb7 details
Key: orb7
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_3sigu_0.1-50.1eV, 1 file(s).
{   'batch': 'ePS co2_2016, batch co2_3sigu_0.1-50.1eV, orbital A2U',
    'event': ' CO2 3sigu-1',
    'orbE': -20.35955918924822,
    'orbLabel': 'orb7',
    'symmetryLabel': '# CO2 3sigu-1'}

*** Job orb6 details
Key: orb6
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_4sigg_0.1-50.1eV, 1 file(s).
{   'batch': 'ePS co2, batch co2_4sigg_0.1-50.1eV, orbital A1G',
    'event': ' CO2 3sigu-1',
    'orbE': -21.709243947049224,
    'orbLabel': 'orb6',
    'symmetryLabel': '# CO2 3sigu-1'}

*** Job stacked details
Key: stacked
No 'job' info set for self.data[stacked].
In [8]:
# # Fix labels for plots
# for key in data.data.keys():
#     data.data[key]['jobNotes']['orbLabel'] = key

System properties¶

Note orbital numbering in table below:

  • Orb is Gamess file output orbital numbering, ordered by energy but not grouped by degeneracy.
  • OrbGrp is grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.
In [9]:
data.molSummary()

*** Molecular structure
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:08:02.095065</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style> 
*** Molecular orbital list (from ePS output file)
EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).
props Sym EH Occ E NOrbGrp OrbGrp GrpDegen NormInt
orb
1 SU -20.6434 2.0 -561.735531 1.0 1.0 1.0 0.974763
2 SG -20.6433 2.0 -561.732810 1.0 2.0 1.0 0.977603
3 SG -11.4445 2.0 -311.420710 1.0 3.0 1.0 1.000000
4 SG -1.5454 2.0 -42.052476 1.0 4.0 1.0 0.999082
5 SU -1.4923 2.0 -40.607552 1.0 5.0 1.0 0.998756
6 SG -0.7978 2.0 -21.709244 1.0 6.0 1.0 0.999632
7 SU -0.7482 2.0 -20.359559 1.0 7.0 1.0 0.999765
8 PU -0.7273 2.0 -19.790841 1.0 8.0 2.0 0.999973
9 PU -0.7273 2.0 -19.790841 2.0 8.0 2.0 0.999973
10 PG -0.5465 2.0 -14.871023 1.0 9.0 2.0 0.999965
11 PG -0.5465 2.0 -14.871023 2.0 9.0 2.0 0.999965 
*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01):
[[1.         0.97476329]
 [2.         0.97760302]]

Compute MF-$\beta_{L,M}$ as function of polarization geometry¶

  • Compute MF results for range of polarization angles $\theta$ (listed as "labels" in plots), set here for a range from $\theta=0$ (z-pol) to $\theta=\pi/2$ (x/y pol).
  • Plot betas vs. Eke and $\theta$.
  • Plot MF-PADs for various channels, single Eke and range of $\theta$. Note that orb assignments may not be correct!

TODO:

  • update code with E-pol calass.
  • confirm assignments and Ekes of relevance.
In [10]:
# With defaults
# eDiag = ep.setPolGeoms(defaultMap='exyDiag')   # Default mapping
# data.MFBLM(RX = eDiag)

# Test alternative z mappings - custom case
tPoints = 10
tRot = np.linspace(0,np.pi/2,tPoints)
pRot = np.zeros(tPoints)
cRot = np.zeros(tPoints)
# labels = tRot.round(2).astype(str)
labels = tRot.round(2)  #.astype(str)
eulerAngs = np.array([labels, pRot, tRot, cRot]).T
eRotRX = ep.setPolGeoms(eulerAngs)

# pRot = [0, np.pi/4]
# tRot = [0, np.pi/2]
# cRot = [0, 0]
# labels = ['x', 'y', 'z']

# Single orb
# orbKeys = 'orb10'
# data.MFBLM(keys=orbKeys, RX = eRotRX)

# All orbs
data.MFBLM(RX = eRotRX)
In [11]:
# Plot all cases with widgets
data.BLMplot(dataType='MFBLM', thres = 1e-2, backend='hv', width=700,  ylim=(-1.5, 2.5), xlim=(0,15))

Applying MFBLM plot defaults, pass overrideMF=True to skip.
BLMplot set data and plots to self.plots['BLMplot']
Out[11]:
True
In [12]:
# Plot single Eke vs. pol rotation ("label" coord)
data.BLMplot(dataType='MFBLM', thres = 1e-2, backend='hv', width=700,  ylim=(-1.5, 2.5),  xDim='Labels', # Erange=[0,2],
            selDims={'Eke':[10.1]})   #, Etype=None)  #,keys='orb10')

Applying MFBLM plot defaults, pass overrideMF=True to skip.
BLMplot set data and plots to self.plots['BLMplot']
Out[12]:
True

X-state (orb10)¶

In [13]:
orbKey = 'orb10'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
            height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels'])   #['Eke','Labels'])

Using default sph betas.
Summing over dims: set()
Plotting from self.data[orb10][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.

Set plot to self.data['orb10']['plots']['MFBLM']['polar']

A-state (orb8)¶

In [14]:
orbKey = 'orb8'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
            height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels'])   #['Eke','Labels'])

Using default sph betas.
Summing over dims: set()
Plotting from self.data[orb8][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.0027680553382225343+5.606882925141018e-17j). This may be unphysical and/or result in plotting issues.

Set plot to self.data['orb8']['plots']['MFBLM']['polar']

B-state (orb7)¶

In [15]:
orbKey = 'orb7'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
            height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels'])   #['Eke','Labels'])

Using default sph betas.
Summing over dims: set()
Plotting from self.data[orb7][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.0002881296903635873+3.586969769800971e-18j). This may be unphysical and/or result in plotting issues.

Set plot to self.data['orb7']['plots']['MFBLM']['polar']

C-state (orb6)¶

In [16]:
orbKey = 'orb6'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
            height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels'])   #['Eke','Labels'])

Using default sph betas.
Summing over dims: set()
Plotting from self.data[orb6][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-1.4300626655083803e-16-1.4065322861916882e-18j). This may be unphysical and/or result in plotting issues.

Set plot to self.data['orb6']['plots']['MFBLM']['polar']

Plot cross-sections and betas¶

These are from ePolyScat's getCro function, and are LF (unaligned ensemble) results. This provides a good, if general, overview.

Overview (all symmetries, length gauge) vs. Eke¶

In [17]:
# NEED TO SET AGAIN AFTER CLASS IMPORT!
import warnings
# warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [18]:
# Comparative plot over datasets (all symmetries only)
Etype = 'Eke'  # Set for Eke or Ehv energy scale
pGauge = 'L'
# Erange=[0, 100]  # Plot range (full range if not passed to function below)
# data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
data.plotGetCroComp(pType='SIGMA', pGauge = pGauge, Etype = Etype)
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:16.354437</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
In [19]:
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype)
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:16.804730</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>

Overview (all symmetries, length gauge) vs. Ehv¶

Same results as above, Ehv scale.

In [20]:
# Comparative plot over datasets (all symmetries only)
Etype = 'Ehv'  # Set for Eke or Ehv energy scale
pGauge = 'L'
# Erange=[0, 100]  # Plot range (full range if not passed to function below)
# data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
data.plotGetCroComp(pType='SIGMA', pGauge = pGauge, Etype = Etype)
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:17.150586</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
In [21]:
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype)
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:17.451630</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
In [22]:
# TODO
# data.plotGetCroComp(pType='BETA', Etype=Etype, backend='hv')

Checks¶

In [23]:
# Betas vs. Gauge
data.plotGetCro(pType='BETA', Etype=Etype)
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:19.109960</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:19.622747</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:20.138269</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"> <cc:Work> <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/> <dc:date>2024-08-29T12:14:20.653255</dc:date> <dc:format>image/svg+xml</dc:format> <dc:creator> <cc:Agent> <dc:title>Matplotlib v3.5.3, https://matplotlib.org/</dc:title> </cc:Agent> </dc:creator> </cc:Work> </rdf:RDF> <style type="text/css">*{stroke-linejoin: round; stroke-linecap: butt}</style>
In [ ]:

Versions¶

In [24]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [25]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
Out[25]:
Thu Aug 29 12:14:21 2024 EDT
OS Linux CPU(s) 64 Machine x86_64 Architecture 64bit
RAM 62.8 GiB Environment Jupyter File system btrfs
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC 11.3.0]
epsproc 1.3.2-dev xarray 2022.3.0 jupyter Version unknown numpy 1.23.5
scipy 1.10.1 IPython 8.13.2 matplotlib 3.5.3 scooby 0.7.2
In [26]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1

* 3d-AFPAD-dev
  dev
  master
93fc0fb4520e1bd89cfebba8f5354b6791c37a1d
In [27]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc

8c9201fdcb98f42875730b94282f50328decab8a	refs/heads/3d-AFPAD-dev
7e4270370d66df44c334675ac487c87d702408da	refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master
baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7	refs/heads/testDev

Env¶

In [28]:
!hostname

8fdb9101787f
In [29]:
!conda env list

# conda environments:
#
base                  *  /opt/conda
In [ ]: